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The effective elastic properties of a unidirectional carbon fiber/epoxy lamina in which 
the carbon fibers are coated with single-walled carbon nanotubes are modeled herein 
through the use of a multi-scale method involving the molecular dy namics/equivalent con- 
tinuum and micromechanics methods. The specific lamina representative volume element 
studied consists of a carbon fiber surrounded by a region of epoxy containing a radially 
varying concentration of carbon nanotubes which is then embedded in the pure epoxy 
matrix. The variable concentration of carbon nanotubes surrounding the carbon fiber 
results in a functionally graded interphase region as the properties of the interphase re- 
gion vary according to the carbon nanotube volume fraction. Molecular dynamics and 
equivalent continuum methods are used to assess the local effective properties of the car- 
bon nanotube/epoxy comprising the interphase region. Micromechanics in the form of 
the Mori- Tanaka method are then applied to obtain the global effective properties of the 
graded interphase region wherein the carbon nanotubes are randomly oriented. Finally, 
the multi-layer composite cylinders micromechanics approach is used to obtain the effec- 
tive lamina properties from the lamina representative volume element. It was found that 
even very small quantities of carbon nanotubes ( 0 . 36 % of lamina by volume) coating the 
surface of the carbon fibers in the lamina can have a significant effect ( 8 % increase) on 
the transverse properties of the lamina (E22, k 23, G 23 and G 12 ) with almost no affect on the 
lamina properties in the fiber direction (£n and 1/12). 

Nomenclature 

A glob Strain concentration factor tensor in global coordinate system 
A loc Strain concentration factor tensor in local coordinate system 
A loc Dilute strain concentration factor tensor in local coordinate system 
Displacement field constant of the j th layer 
C Effective stiffness tensor in global coordinate system 
Cf Fiber volume fraction 

C g } ob Fiber stiffness tensor in global coordinate system, Pa 
Cj Fiber stiffness tensor in local coordinate system, Pa 
C m Matrix stiffness tensor, Pa 
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Cn Axial stiffness component for transversely isotropic material, Pa 
En Axial Young’s modulus for transversely isotropic material, Pa 
E 22 Transverse Young’s modulus for transversely isotropic material, Pa 
I Fourth order identity tensor 

N Total number of layers in the composite cylinder assemblage, m 
L Length of composite cylinder assemblage 

r Radial position in cylindrical coordinates 

rj Radial layer boundary, m 

S loc Eshelby tensor in local coordinate system 
uP Radial displacement component of the j th layer, m 
uP z-direction displacement component of the j th layer, m 
^-direction displacement component of the j th layer, m 
V Volume, m 2 3 

V 0 Volume fraction of carbon nanotubes in Region II 

Vo Initial volume fraction of carbon nanotubes in Region III 

W Volume averaged strain energy of the composite cylinder assemblage, N/m 4 

IF ( eff ) Volume averaged strain energy of the effective material, N/m 4 

Xi Global coordinate direction 

Xi Local coordinate direction 

A Thickness of Region III, m 

£0 Applied strain 

ep Small strain components of the j th layer 

k >23 Plane strain bulk modulus for transversely isotropic material, Pa 
X i Lame constant for isotropic phase, Pa 

fii Lame shear constant for isotropic phase, Pa 
fii 2 Axial shear modulus for transversely isotropic material, Pa 
fi 23 Transverse shear modulus for transversely isotropic material, Pa 
V 12 Axial Poisson’s ratio for transversely isotropic material 
Cauchy stress components of the j th layer, Pa 

Subscripts and Superscripts 

1 Variable number 

j Variable number 

(j) Variable layer designation in composite cylinders model 

k Variable number 

r Radial component in cylindrical coordinates 

2 z-direction component in cylindrical coordinates 

9 ^-direction component in cylindrical coordinates 


I. Introduction 

A pplications for carbon nanotubes (CNTs) in aerospace vehicles are likely to focus in the near term 
on the enhancement of the mechanical, thermal and electrical properties of materials currently in use. 
The exceptional properties of CNTs and their relatively low density are believed to be able to provide 
large improvements in properties with only small quantities of CNTs. For example, enhancements in elastic 
Young’s modulus of unfunctionalized CNTs in CNT-epoxy composites of up 20% in tension and 75% in 
compression have been observed for CNT volume fractions of just 5%. 1 For functionalized CNTs, that 
is CNTs which have been chemically altered to better interact with a given matrix, increases in Young’s 
modulus of 30% and 60% have been observed at volume fractions of 1 and 4%, respectively. 2,3 As such, one 
proposed implementation of CNTs is as structural enhancement of traditional carbon fiber/epoxy composite 
laminates, as shown schematically in Figure 1, where it is believed that the selective use of CNTs as surface 
treatments can improve interface strength between the carbon fibers and epoxy, and thereby improve the 
fracture toughness of the composite. 
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Randomly Oriented CNTs 
in Regions II and III 



Molecular Dynamics 
CNT-Epoxy RVE 


(c) 


(d) 


Figure 1. Motivation: CNT-Sprayed Carbon Fiber Composite Analysis Scales (a) Lamina Scale (mm); (b) 
Lamina RVE Scale (//m); (c) Region II and III Scale ( fim *); (d) CNT RVE Scale (nm) 


As a result of the inclusion of nanoscale materials, analysis of such composites, often termed nanocom- 
posites, requires the development of multi-scale models which can span the inherent nanometer to millimeter 
and beyond length scales. For example, in the composite laminate schematic of Figure 1, the carbon fiber 
lamina has dimensions on the order of millimeters, the lamina representative volume element (RVE) (Fig- 
ure 1(b) is on the order of microns, the CNT coating surrounding a carbon fiber (Figure 1(c)) on the order of 
a tenth of a micron, and finally the CNT RVE in Figure 1(d) is on the order of nanometers. Modeling across 
such disparate length scales necessitates the use of tools such as molecular dynamics (MD), micromechanics, 
composite laminate analysis, and finite element modeling. 

Multi-scale models have been applied using both atomistic simulations and micromechanics to assess the 
constitutive properties of various functionalized nanotube materials. 4-9 Other researchers have used atomistic 
simulations of epoxy /nanotube composites in order to address nanotube pull-out from the epoxy. 10 A more 
recent study on nanotubes chemically bonded into the epoxy matrix showed a Young’s modulus of up to 160 
GPa in the direction of the nanotube axis, and 4-8 GPa in the transverse direction at a nanotube volume 
fraction of 25%. 11 

In the present work, the constitutive properties of a unidirectional carbon fiber/epoxy lamina in which 
each carbon fiber is coated with CNTs are investigated. The large difference in scales between the carbon 
fiber, with radius on the order of microns, and the carbon nanotubes, with radius on the order of nanometers, 
suggests the use of a multi-scale approach. The multi-scale approach herein makes use of molecular dynamics 
simulations and the equivalent continuum method on the CNT scale and continuum micromechanics analysis 
on the carbon fiber interphase and lamina scales in assessing effective lamina elastic properties. The objective 
is to calculate the effective elastic constants of unidirectionally-reinforced carbon fiber/epoxy laminae, where 
the CNT distribution around the fiber is varied from a high concentration near the carbon fiber surface 
to zero CNT concentration as a result of the spraying the process. The constitutive properties of the 
nanoscale components are modeled with the molecular dynamics-equivalent continuum model. The CNT 
orientation and distribution in volume fraction are modeled using the results of the molecular dynamics- 
equivalent continuum model in a Mori-Tanaka micromechanics approach to calculate the effective properties 
of randomly oriented CNTs in epoxy. Finally, the multi-layer composite cylinders method is applied to the 
lamina representative volume element to obtain the effective lamina properties for carbon fibers in epoxy 
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surrounded by a graded interphase of varying CNT concentration. A description of each stage of the multi- 
scale model is presented, and the results of applying the multi-scale approach are given for several types of 
CNT coatings in the lamina. 


II. Micromechanics Modeling 

A. Multi-layer Composite Cylinders Model 

The composite cylinders method was originally proposed by Hashin and Rosen 12 to obtain the effective elastic 
constants of unidirectional fiber reinforced composites. In this method, a set of five elasticity solutions are 
obtained for both a composite cylinder assemblage (cylindrical fiber surrounded by cylindrical matrix) and 
for a transversely isotropic effective homogeneous cylinder. The strain energy of the composite cylinder 
assemblage is equated to that of the effective cylinder, and the five independent effective elastic constants 
are determined. The composite cylinders method of Hashin and Rosen worked well for four of the five 
independent elastic constants of the effective transversely isotropic material, but could only provide bounds 
on the fifth, the effective transverse shear modulus 
Christensen and Lo 13 proposed an alternate 
model for the effective transverse shear modulus in 
the form of the generalized self-consistent compos- 
ite cylinder in which the original composite cylinder 
assemblage is surrounded by an additional concen- 
tric cylinder having the same elastic properties as 
the effective homogeneous cylinder. The generalized 
self-consistent composite cylinder solution resulted 
in an explicit solution for the effective transverse 
shear modulus that fit within the Hashin and Rosen 
bounds. The methods of Hashin and Rosen and 
Christensen and Lo were later consolidated and and 
generalized in applications involving coated fiber 
composites 14, 15 resulting in a multi-layer composite 
cylinders method (Figure 2). Recently, 16 multi-layer 
composite cylinders models have been applied to 
carbon nanotube reinforced epoxies wherein carbon 
nanotubes have been treated as continuum sheets of 
rolled graphene surrounded by interphase layers to 
simulate differing amounts of load transfer as a re- 
sult of functionalization, but have lacked direct MD 
coupling. 

In the present work, the effective properties for 
the carbon fiber lamina RVE are obtained using the 
multi-layer composite cylinder. As shown in Fig- 
ure 1(b), the layers are the carbon fiber (Region I), the high concentration of carbon nanotubes interphase 
region (Region II), the variable concentration of carbon nanotubes interphase region (Region III), and the 
epoxy matrix (Region IV). As an example, the procedure for determining the effective plane strain bulk 
modulus (^ 23 ) for the lamina RVE is provided. Each layer of the composite cylinder assemblage is taken to 
have the following displacement field (in cylindrical coordinates with x\ — z): 

1 4/> = B[ j) r + B ( 2 j) - 

w]/ 1 = 0 for I'j _ 1 <r <rj (1) 

— 0 

where the superscript (]) denotes the region (either I, II, III or IV) and b[^ and are constants. The 
range of each layer is defined by the inner and outer radii, rj- ± and r^, respectively, with ro = 0 for 
Region I and 77 y for Region IV dependent on the carbon fiber volume fraction. The strains in each layer 
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Composite Cylinder 
Assemblage 


A 


Homogeneous 
Effective Material 


Figure 2. Schematic representation of the composite 
cylinder method. The composite cylinder assemblage 
(left) consists of N concentric cylinders whose strain 
energy under load is equivalent to that of an effective 
homogeneous material (right). 
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are determined from the displacement fields using the small strain strain-displacement relations, with the 
layer stresses determined from the appropriate constitutive relation for each layer (isotropic linear elastic for 
Regions I-IV). The constants b[^ and are determined through the application of boundary conditions 
and interface matching conditions. The boundary conditions for determining the plane strain bulk modulus 
are given by: 

= S 0 r N (2) 

Bi 1] = o (3) 

where is the radial strain applied at the boundary and Eqn. (3) insures that the displacement in Region 
I is bounded at r$. The matching conditions at the layer interfaces impose the continuity of displacement 
and traction conditions and are again given by: 


,0)| =7/( /e + 1 )| 


r(*)| 


= rr( k + 1 ) I 


rr \r=r k — ^ rr \r=r k 


( 4 ) 

( 5 ) 


where k ranges from one to N — 1. Thus, we have a system of 2N equations to solve for the 2N unknown 
constants in the displace field, b\^ . 

The elasticity solution for the effective homogeneous cylinder follows the same procedure. The displace- 
ment field in the homogeneous cylinder is taken to be given by: 


4 eff) = B{ eS) r + B ( 2 eS) ^ 

Ug eS ^ = 0 for 0 < r < rjv ^ 


,(eff) _ 


= 0 


where the superscript (eff) denotes the effective homogeneous material. Again the strains are determined 
from the strain-displacement relations and the stresses from the constitutive relations which, for the effective 
homogeneous cylinder, are the unknown effective constants of a transversely isotropic material. The boundary 
conditions are the same as those applied in Eqns. (2) and (3), and no matching conditions are necessary as 
the material is homogeneous. 

Having solved for all of the unknown constants in the displacement field, the volume averaged strain 
energies of the composite cylinder assemblage, W, and of the effective homogeneous cylinder, W^\ can be 
obtained from: 

W — (o’ijSij) (7) 

and 

^ (eff) = <df Mf *> ( 8 ) 

where in both equations the summation convention is invoked, and where the lack of superscripts in Eqn. (7) 
is intended to indicate that the stress and the strain in the brackets is for the composite cylinder assemblage 
as a whole, i.e, across all layers. The brackets denote the volume average as follows: 

(•) = y I • dV (9) 

V 

It can be shown that volume averaged strain energy of the effective homogeneous cylinder can be expressed 
as: 

W {eS) = 4 K 23 e 2 0 (10) 

The volume averaged strain energy of the composite cylinder assemblage can be expressed in terms of the 
stress and strain in the layers as: 

w= v\hUU < u > 
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The effective plane strain bulk modulus is then determined from equating the volume averaged strain 
energies of the composite cylinder assemblage and the effective homogeneous cylinder, i.e.: 


w = W eS 


( 12 ) 


It can be shown that Eqn. ( 12 ) can be solved for ^23 which can be expressed in the following form which 
depends only on the stress and displacement of the N th phase: 


ft 23 


7 W, 

' rr | r=rN 


( 'ui N) \ r=rN /r N ) 


( 13 ) 


which is analogous to the definition of the plane strain bulk modulus given by ^23 = ^22/2622- Eqn. ( 13 ) 
can then be expressed in terms of the now known constants as: 


ft 23 


(/in + Atv )r 2 N B^ — unB^ 


A N ) 


• N 


b[ N) + B. 


(TV) 


( 14 ) 


where it is has been assumed that the N th phase is isotropic (fijsr and Xjy are the Lame constants of the N th 
phase, i.e., Region IV). Similar procedures are followed for obtaining the remaining four transverse isotropy 
constants: the axial Young’s modulus, En, the axial stiffness component, Cn, the axial shear modulus, /ii2, 
and the transverse shear modulus, /i23 - 15,16 From these five, other properties such as the transverse Young’s 
modulus, E 2 2, and the axial Poisson’s ratio, v 12 can be calculated. 


B. Effective Properties of Randomly Oriented Fibers 

In order to identify constitutive properties for Re- 
gions II and III of the lamina RVE in Figure 1 (b), 
it is necessary to consider the effective properties of 
randomly oriented CNTs. Current MD simulations 
would find it computationally difficult to model an 
RVE of sufficient size to consider all possible ori- 
entations of CNTs needed to capture the effects of 
random orientation. As such, micromechanics tools 
in the form of the Mori- Tanaka method 17, 18 are em- 
ployed in a randomization procedure. 

The randomization procedure is shown schemat- 
ically in Figure 3 . Inclusions, in this case CNTs, at 
each distinct orientation relative to the global coor- 
dinate system {X\) of the composite are considered 
as separate phases in a multi-phase Mori- Tanaka ap- 
proach. However, in the local inclusion coordinate 
system (a^), the inclusions are identical and there- 
fore have the same strain concentration factor in the local coordinate system given by: 



of Inclusions 



Local Inclusion 
Orientation 


Figure 3. Schematic representation of the randomiza- 
tion procedure. Each inclusion’s stiffness and concen- 
tration factor in the local coordinate system (x{) are 
rotated to the global coordinate system (X{) and aver- 
aged as a separate phase. 


A loc = A 


loc 


(1 - Cf)I + CfA 


loc 


-1 


( 15 ) 


where A loc is the Mori- Tanaka dilute strain concentration factor which can be expressed as: 

A loc = [l + S loc (C TO )- 1 (C 1 / oc -C ra )] _1 (16) 

where Cf is the phase volume fraction, I is the fourth order identity tensor, S is the Eshelby tensor 19 (in this 
case for a circular cylinder), and C m and Cf are the matrix and fiber stiffness tensors, respectively. The 
superscript loc in Eqns. (15) and (16) indicates that the given quantity is expressed in the local coordinate 
system, and it is noted that the matrix is assumed to be isotropic so that the local and global coordinate 
representations of the matrix stiffness (C m ) are the same. 
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The orientation of the local coordinate system of a given phase relative to the global coordinate system 
of the composite is expressed in terms of general Euler angle rotations. Subsequently, a general change 
of basis is applied to rotate the local stiffness tensor and concentration factor into their global coordinate 
system representations, C® lob and A glob , respectively. The effective properties for the composite in the global 
coordinate system can then be expressed as: 

C = C m + c/ (cf b - C m )A glob (17) 

Finally, the effective properties of the composite as a whole are then obtained by integrating the second term 
of Equation 17 over all orientations, and are found to isotropic. 

III. Application: CNT-sprayed Carbon Fiber Lamina 

The multi-layered composite cylinder method 
and the micromechanical randomization procedures 
described above are applied in the modeling of CNT- 
sprayed unidirectional carbon fiber-epoxy compos- 
ites (Figure 1). The lamina, Figure 1(a), is con- 
structed such that the carbon fibers are first coated 
with CNTs, and then cured in epoxy. The schematic 
in Figure 1(a) also identifies the principal material 
axes with respect to the fiber direction, with X\ 
referred to as the longitudinal, or axial direction, 
and X 2 and V3 referred to as the transverse, or in- 
plane directions. The epoxy of interest is comprised 
of bisphenol F resin reacted with epichlorohydrin 
and cured with triethylenetetramine (TETA). The 
CNTs sprayed on the carbon fiber were taken as 
(10,10) single-walled carbon nanotubes with a ra- 
dius 0.678 nm. As the CNTs are not functionalized, 
nor are they grown on the carbon fiber surface, it 
is assumed in the model that there is no chemical 
bonding between the 3 components: carbon fiber, carbon nanotubes, and the epoxy. 

The effective properties of the lamina in Figure 1(a) are determined through the application of the 
multi-layer composite cylinders method on the lamina RVE shown in Figure 1(b) (Xl - axis, positive out- 
of-the-page) . The RVE identifies four regions. Region I is the carbon fiber itself, where rp represents the 
carbon fiber radius (taken as 5/ira). Region II represents a region of relatively high CNT volume fraction 
very near the carbon fiber surface, and ranges from the carbon fiber radius at rp to ru for a thickness of 
0.004 microns. Region III represents a much larger CNT rich region surrounding the carbon fiber where the 
CNT volume fraction decays from a value close to that of Region II to a zero volume fraction of CNTs. Both 
Region II and III are comprised of randomly oriented CNTs and epoxy (Figure 1(c)). The width of Region 
III ranges from ru to rm, and is denoted by the parameter A, where A = rm — ru and takes a value of 
0.3 fim. Region IV is the remaining epoxy matrix (i.e., no CNTs), and ranges from rm to a maximum radius 
denoted by the parameter r\\ which is determined based on the carbon fiber volume fraction. 

A detailed graphical representation of the CNT volume fraction distribution versus radial distance from 
the carbon fiber center for each region of the lamina RVE is shown in Figure 4. Regions I and IV, the carbon 
fiber and the epoxy matrix, respectively, contain no CNTs and thus have a CNT volume fraction of zero. 
The CNT volume fraction in Region II is constant and is defined by the parameter V 07 shown in Figure 4, 
which is taken to be 15.55%. The CNT volume fraction in Region III decays from a value Vo to zero as the 
radial distance from the carbon fiber increases. The method of decay for the CNT volume fraction in Region 
III is given by Eqn. (18): 

v f (r) = v 0 e- a{r ~ r ^ - voe~ aA (18) 

where a — 9, v 0 — 0.1535, and A = 0.3 fim is radial thickness of Region III. It should be noted that the 
thicknesses and distribution of CNT volume fractions in Regions II and III are assumed values for modeling 



Radial distance from carbon fiber, r(n) 


Figure 4. Distribution of CNT volume fraction around 
the carbon fiber for the Lamina RVE in Figure 1(b). 
Region II represents a region of higher CNT volume 
fraction near the carbon fiber surface. Region III rep- 
resents a larger region of decaying CNT volume fraction 
which decays according to Eq. (18). 
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purposes as at present there is insufficient characterization of the fabricated composites to determine these 
values experimentally. 

With the microstructure of the lamina RVE identified, the procedure for determining the effective lamina 
properties is outlined as follows. First, MD simulations are used to generate representative volume elements 
of the of the CNT-epoxy constituents of Regions II and III (Figure 1(d)). Second, the equivalent-continuum 
model 20 is used in conjunction with these MD simulations to calculate the continuum level elastic constants 
of the CNT-epoxy RVEs. Next, the elastic constants of randomly oriented CNTs in epoxy are calculated 
as a function of CNT volume fraction using the Mori- Tanaka method. The elastic constants determined 
from Mori- Tanaka randomization are then mapped onto the CNT volume fraction distribution of Figure 4 to 
obtain the distribution of elastic constants as a function of radial distance from the carbon fiber for Regions II 
and III. Finally, the multi-layer composite cylinders method is applied to the lamina RVE using the mapped 
distribution in elastic constants. These steps are discussed in detail below, as are results obtained from 
the applying them towards determining the effective properties of the CNT-sprayed carbon fiber composite 
lamina described above. 

A. Molecular Dynamics-Equivalent Continuum CNT-Epoxy Modeling 

Molecular dynamics simulations of CNTs 
in epoxy are performed for the molecu- 
lar RVE shown in Figure 1(d). The RVE 
of the molecular structure of the mate- 
rial system was generated via initial en- 
ergy minimization in the MD simulations 
and is reflective of a CNT volume frac- 
tion of approximately 15%. The molecu- 
lar structure of the epoxy was a network 
of bisphenol F resin molecules reacted 
with epichlorohydrin and cross-linked via 
triethylenetetramine (TETA) molecules. 

The TETA molecule is able to cross-link 
with up to four resin chains, and on av- 
erage 2.3 of the 4 TETA amine groups in 
the resultant structure were involved in 
cross-linking. The equilibrium structure 
of the CNT-epoxy RVE was obtained by 
compressing the epoxy network until the minimum potential energy configuration at finite temperature (300 
K) was reached. 

The CNT-epoxy RVE was simulated with the AMBER force field 21 using the LAMMPS 22 and DL- 
POLY 23 simulation codes. All of the CNT and graphite parameters were taken from the AMBER force 
field. Several of the force field parameters of the epoxy network were derived from ab initio calculations 
carried out with the NWChem package. 24 The partial atomic charges were derived from an RESP fit to 
a Hartree-Fock calculation of wavefunction. The bond lengths, angles, and dihedrals, were calculated by 
optimizing the geometry of parts of the epoxy network at the HF/STO-3G level. Partial charges were 
then fit at this geometry from a HF/6-31G* calculation. The rest of the parameters, were taken from the 
AMBER force field. In the simulations, electrostatic interactions were calculated with Coulomb’s law, and 
an Ewald summation was used for the long range Coulombic interactions. Electrostatic interactions and 
1-4 non-bonded dihedral interactions were scaled by the standard AMBER parameters of 0.83 and 0.50 
respectively. 

In the MD simulations the molecular structure RVE is subjected to periodic boundary conditions in 
all three directions. The displacement fields from Ref., 7 which are consistent with those used in obtaining 
continuum level elastic properties, were applied to the RVE structure in strain increments of 0.0025 every 
10,000 steps at 1 fs each (10 ps per increment) up to a maximum strain of greater than 0.02. The application 
of the displacement field included displacing both the periodic boundaries and the atomic positions, then 
equilibrating the system with MD to its new state point. The elastic constants were then obtained for the 
molecular RVE from equating the energies of deformation for the molecular system, which were calculated 
from molecular dynamics simulation, to the energies of deformation of an equivalent continuum homogeneous 
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Table 1. The elastic constants ( GPa ) obtained from the MD- 
equivalent continuum method and the corresponding engineering 
constants ( GPa , except V 12 ) obtained from Eqns. (19) through 
(24). 


MD-Equivalent Continuum 

Graphite CNT 

Elastic Constants 

Engineering Constants 

C n 

99.00 

En 

93.44 

109.30 

C22 

12.20 

E22 

0.40 

2.28 

C 12 

8.20 

V \2 

0.34 

0.37 

£ 

CO 

12.00 

K 23 

12.10 

3.83 

C44 

0.20 

G 23 

0.10 

0.67 

c 66 

0.33 

Gl2 

0.16 

0.69 
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solid under identical deformations. 

The effective stiffness tensor components for the CNT-epoxy system obtained using the MD /equivalent 
continuum methods are provided in Table 1. Also provided in Table 1 are the expression of the stiffness 
tensor in terms of the engineering elastic constants for transversely isotropic materials, with the relationships 
between the two representations given in Eqns. (19) through (24). 


En = C n - 


C h 

C22 + C23 


(19) 


, , Cf 2 (2C 23 - C22) - ClgCll 

22 T /o /o si 2 

^11^22 — ^12 

(20) 

Cn 

C22 + C23 

(21) 

ft 23 = 2 (^22 + C23) 

(22) 

G 23 = \ (C 22 - C 23 ) 

(23) 

Gn = \c 66 

(24) 



Composite Cylinder 
Assemblage for CNT-Epoxy 


Figure 5. Composite cylinder assem- 
blage used in initial micromechanics 
evaluation of the effective properties 
of CNT-epoxy composites. 16 The re- 
sulting effective properties of this com- 
posite cylinder assemblage are com- 
pared with those obtained from MD- 
equivalent continuum correspondence 
of the MD-RVE in Figure 1(d) in Ta- 
ble 1. 


For comparison purposes, the effective properties of CNTs em- 
bedded in epoxy at 15% volume fraction obtained using the com- 
posite cylinders solution are also provided in Table l. 16 In these so- 
lutions, the composite cylinder assemblage consists of a hollow and 
solid portion of the CNT surrounded by the epoxy matrix as shown 
in Figure 5. Here the composite cylinder method describe above is 
applied for N — 2, with the hollow region accounted for through the 
application of traction free conditions on the CNT inner surface. As 
such, the only difference for the composite cylinders method is the 
change of the boundary condition in Eqn. (3) to al r \ r=ro — 0, with 
ro equal to the CNT inner radius of 0.54nra . The CNTs are thus 
represented as a rolled sheet of graphene, with the solid portion of 
the CNT having isotropic properties reflective of the in-plane prop- 
erties of graphene (E — 1100 GPa 7 v — 0.14). 25 Table 1 indicates 
that the composite cylinders evaluation of the effective CNT-epoxy 
properties differs greatly from the MD-equivalent continuum values, 
especially for E 2 2 7 ft23, G \2 and G23. It is of interest to note that 
these four properties have been observed to be largely affected by the 
presence of interphase regions. 16 The MD simulations used to obtain 
the effective CNT-epoxy properties would automatically and accu- 
rately include the presence of an interphase region, which is believed 
to be the source of the large differences observed in Table 1. 

B. Random Orientation of Effective CNT-Epoxy Proper- 
ties 


The CNTs in Regions II and III are assumed to be randomly ori- 
ented. As such, the effective CNT-epoxy properties from the MD- 
equivalent continuum method provided in Table 1 are used to define 
the local stiffness and Eshelby tensor of Eqn. (16). Here the Eshelby tensor applied is for a circular cylinder 
in an isotropic matrix as the surrounding epoxy is considered to be isotropic linear elastic. The global aver- 
aging procedure of Eqn. (17) was applied to obtain the stiffness tensor for the randomly oriented CNT-epoxy 
system in Regions II and III. Despite being locally transversely isotropic, the effective properties obtained 
from the global averaging process are isotropic. As such, the effective Young’s modulus (E) and Poisson’s 
ratio (y) for Regions II and III as a function of CNT volume fraction are obtained and provided in Fig- 
ures 6(a) and 6(b), respectively. The Young’s modulus and Poisson’s ratio versus CNT volume fraction are 
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then mapped into the CNT volume fraction distribution of Figure 4 so that the material properties and 
material property variation of Regions II and III are identified. 


Randomly Oriented Effective Poisson's Ratio 

0.45 
0.43 
0.41 
0.39 
0.37 
0.35 

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 

Volume Fraction of CNTs 

(bF 

Figure 6. Effective elastic properties for randomly oriented MD-equivalent continuum CNT representations 
obtained from Eq. (17) and used for Regions II and III of Figure 1(b) corresponding to the CNT volume 
faction of Figure 4. (a) Effective Young’s modulus; (b) Effective Poisson’s ratio. 




(a) 


C. Composite Cylinders Solution for Effective Lamina Properties 

With the material properties of Regions II and III identified, the multi-layer composite cylinders method can 
be applied to determine the effective properties of the carbon fiber lamina. The effects of the sprayed CNTs 
on the effective lamina properties are studied in five cases below. The first two are single interphase cases, 
and the remaining three are multiple interphase cases. Detailed descriptions of each case and discussions on 
the resulting effective properties are provided below. 

1 . Single Interphase Solutions 

Two single interphase cases, that is a composite cylinder assemblage of N — 3, are considered in an effort to 
demonstrate the different effects of Regions II and III. In the first case, the effects of having the carbon fiber 
in the lamina treated with a small region with a relatively high volume fraction of CNTs are investigated. 
This case is representative of having Region I, Region II, and Region IV only. In the second case, the 
effects of having the carbon fiber in the lamina treated with a larger region of CNTs, but with a variable 
volume fraction of CNTs are is investigated. However, in this case, which is representative of having Region 
I, Region III, and Region IV only, the properties for Region III are obtained from the properties of the 
random orientation of effective CNTs at the CNT volume fraction corresponding to the average volume 
fraction of Region III in Figure 4 (i.e. 4.86%). The material property and geometry data for both cases 
are summarized in Table 2. It should be noted that the epoxy elastic constants provided for Region IV in 
Table 2 were obtained from MD-equivalent continuum simulations similar to those done on the molecular 
RVE in obtaining the effective CNT-epoxy properties; the only difference being the lack of a CNT in the 
molecular RVE. Typical values for storage moduli of the neat cured epoxy are 2-3 GPa at 250-300 K. 26 

In Figure 7(a), the effect of the single interphase regions in Cases 1 and 2 on the effective axial Young’s 
modulus (En) of the lamina are compared to the axial Young’s modulus results for a carbon fiber lamina 
with no CNTs (i.e., a composite cylinder assemblage of N — 2 where the two phases are Region I and Region 
IV). From Figure 7(a) it is observed that the effect of the interphase region on the effective axial modulus 
is negligible. For the no CNT case, the axial modulus is known to be well approximated by the rule of 
mixtures. The axial Young’s modulus for Cases 1 and 2 are similarly well approximated by a three-phase 
rule of mixtures, however, given the large differences between the fiber and matrix properties, and that the 
interphase regions are of the same order as the matrix properties, it is difficult to observe any difference 
between the no CNTs and single interphase cases. In fact, at carbon fiber volume fractions typically found 
in lamina, e.g., 60% carbon fiber, the percent difference in axial Young’s modulus between Case 1 and the no 
CNT case is 0.0017%. For Case 2, the percent difference is 0.041%. Although not shown here, the effective 
Poisson’s ratio of the lamina is also marginally effected by the interphase regions in Cases 1 and 2. Unlike 
the axial modulus, the Poisson’s ratio consists of both axial and transverse strain components resulting in 
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Table 2. The Young’s modulus and Poisson’s ratio of each phase used for two single interphase composite 
cylinder solutions for CNT-sprayed carbon fiber lamina. The first case is indicative of a small region of a 
relatively high concentration of CNTs near the carbon fiber. The second case is indicative of a larger lower 
CNT concentration region representing the average concentration of Region III CNT volume fraction decay 
of Figure 4. 



Region 

Radial Distance jim 

Young’s Modulus GPa 

Poisson’s Ratio 

Case 1 

Region I 

0-5.0 

230.00 

0.3000 

Region II 

5.0 - 5.004 

3.8608 

0.3696 

Region IV 

> 5.004 

1.4400 

0.4200 

Case 2 

Region I 

0-5.0 

230.00 

0.300 

Region III 

5.0 - 5.3 

2.2039 

0.4003 

Region IV 

> 5.3 

1.4400 

0.4200 


slightly higher percent differences for Cases 1 and 2 relative to the no CNT case, e.g., 0.015% and 0.47%, 
respectively, at 60% carbon fiber volume fraction. 



Figure 7. Effective Axial Young’s modulus (En) (a) and plane strain bulk modulus (^ 23 ) (b) of CNT-sprayed 
carbon fiber lamina. The label No CNTs refers to a composite cylinders assemblage for the Lamina RVE con- 
sisting of just Regions I and IV. The labels CNTs Near C-Fiber and CNTs Far from C-Fiber refer single interphase 
cases where the single interphase corresponds to Region II and Region III, respectively, with the Region 
III interphase value taken from an average of the CNT volume fraction decay in Figure 4. The number in 
parenthesis denotes number interphase regions. 


The effective plane strain bulk modulus (^23) results for Cases 1 and 2 are compared to the no CNT case 
in Figure 7(b). It is observed that, despite having a higher volume fraction of CNTs in the single interphase 
region of Case 1, the small thickness of the interphase renders the contribution of the CNTs to the effective 
plane strain bulk modulus to be negligible (0.11% difference relative to the no CNT case). In contrast, 
the low CNT volume fraction but large thickness of the interphase region in Case 2 produces observable 
differences relative to the no CNT case, particularly at high carbon fiber volume fractions. At 60% carbon 
fiber volume fraction, the percent difference in plane strain bulk modulus for Case 2 relative to the no CNT 
case is 4.32% and becomes as large as 11.5% at 80% carbon fiber volume fraction. Reasons for the differences 
between Cases 1 and 2 can be explained in terms of the total lamina CNT volume fraction. The high CNT 
volume fraction, small interphase thickness of Case 1 results in a total lamina CNT volume fraction of 0.015% 
at 60% carbon fiber volume fraction. a In contrast, the low CNT volume fraction but large thickness of the 
interphase in Case 2 results in a lamina CNT volume fraction of 0.36% at 60% carbon fiber volume fraction. 
Though not shown here, the remaining effective lamina properties, G12, G 23, and ^22 > show very similar 

a Note that the properties and thicknesses of the single interphase regions remain fixed as carbon fiber volume fraction varies. 
Increasing the carbon fiber volume fraction reduces the epoxy (Region IV) volume fraction by making r iv smaller until it equals 
the interphase radius. As such, at higher carbon fiber volume fractions, the interphase lamina volume fraction is increased and 
therefore, so to is the lamina CNT volume fraction. 
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trends to the plane strain bulk modulus results provided in Figure 7(b). The only difference between these 
properties and the plane strain bulk modulus being in the magnitude of percent difference relative to the no 
CNT case which remains at less than 0.3% for Case 1 and is on the order of 7-8% for Case 2 at 60% carbon 
fiber volume fraction. As such, for the multiple interphase results which follow, only the plane strain bulk 
modulus will be provided. 

2. Multiple Interphase Solutions 


Table 3. The Young’s modulus and Poisson’s ratio of each phase used for three multiple interphase composite 
cylinder solutions for CNT-sprayed carbon fiber lamina. These cases are indicative of approximations of the 
CNT concentration region of Regions II and III CNT volume fraction decay of Figure 4. 



Region 

Radial Distance \im 

Young’s Modulus GPa 

Poisson’s Ratio 

2-Interphase 

Region I 

0 - 5.0 

230.00 

0.3000 

Region II 

5.0 - 5.004 

3.8608 

0.3696 

Region III 

5.004 - 5.304 

2.2039 

0.4003 

Region IV 

> 5.304 

1.4400 

0.4200 

3-Interphase 

Region I 

0 - 5.0 

230.00 

0.3000 

Region II 

5.0 - 5.004 

3.8608 

0.3696 

Region IH.a 

5.004 - 5.154 

2.7286 

0.3880 

Region IH.b 

5.154 - 5.304 

1.6654 

0.4128 

Region IV 

> 5.304 

1.4400 

0.4200 

4-Interphase 

Region I 

0 - 5.0 

230.00 

0.3000 

Region II 

5.0 - 5.004 

3.8608 

0.3696 

Region III. a 

5.004 - 5.104 

3.0036 

0.3826 

Region IH.b 

5.104 - 5.204 

2.0096 

0.4030 

Region III.c 

5.204 - 5.304 

1.5700 

0.4157 

Region IV 

> 5.304 

1.4400 

0.4200 


Such large differences between Cases 1 and 2 in- 
dicate a need to accurately represent the interphase 
region in determining the effective lamina proper- 
ties. As such, more accurate representations of the 
graded interphase region of CNT treated carbon 
fibers are also investigated through the use of mul- 
tiple interphase regions. In these cases, all four re- 
gions are considered, i.e., Region I, Region II, Re- 
gion III, and Region IV. The properties for Region 
II are obtained as previously discussed for Case 1 
above. The graded interphase of Region III is repre- 
sented using increasing numbers or piecewise contin- 
uous subregions. The subregion properties are ob- 
tained by taking the average Young’s modulus and 
Poisson’s ratio over the subregion thickness from the 
mapping of the effective properties of randomly ori- 
ented CNT-epoxy RVE into the CNT volume frac- 
tion distribution of Figure 4. The material property 
and geometry data for each multiple interphase case 
are summarized in Table 3, with the 4-Interphase 
Case distribution in Young’s modulus shown graph- 
ically in Figure 8 and compared with the exponential decay of Eqn. (18). These three multiple interphase 



Figure 8. Radial distribution in Young’s modulus in Re- 
gions II and III obtained from the correlation between 
the desired CNT volume fraction distribution given in 
Figure 4 and the random orientation Young’s modulus 
results of Figure 6(a). Region III is approximated by 
three piecewise continuous subregions which, with Re- 
gion II, results in a composite cylinder assemblage of 
N = 6. 
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cases are represented by composite cylinder assemblages of N — 4, 5, and 6. The 2-Interphase case represents 
a combination of the two single interphase cases previously discussed. The 3-Interphase and 4-Interphase 
cases are intended to demonstrate increasingly accurate representations of the graded interphase of Region 

III. 

Effective plane strain bulk modulus results at 
high carbon fiber volume fractions (greater than 
80%) are provided for all cases in Figure 9. The re- 
sults in Figure 9 show that the no CNT case brackets 
the other cases from below while the 2-Interphase 
case brackets the results from above. The order 
of the resulting plane strain bulk moduli from the 
bottom up is No CNTs, Case 1, 4-Interphase, 3- 
Interphase, Case 2, and 2-Interphase, which at 80% 
carbon fiber volume fraction amounts to percent dif- 
ferences relative to the No CNT case of 0%, 0.28%, 

9.85%, 10.19%, 11.49% and 11.85%, respectively. It 
is of interest to point out that the order of increas- 
ing total lamina CNT volume fraction does not ex- 
actly correspond to the order of increasing plane 
strain bulk modulus. The total lamina CNT vol- 
ume fractions are 0%, 0.020%, 0.47%, 0.49%, 0.48% 
and 0.50% for the No CNTs, Case 1, 4-Interphase, 

3-Interphase, Case 2, and 2-Interphase cases, re- 
spectively. The 3-Interphase case actually has a 
higher total lamina CNT volume fraction than Case 
2 which has the higher plane strain bulk modulus of 
the two. This indicates the importance of the dis- 
tribution of the CNTs around the carbon fiber. In 
Case 2 there is a uniform distribution of the CNTs 
in the interphase region, but in the 3-Interphase the 
distribution of CNTs into three separate subregions 
results in an average Young’s modulus for the Re- 
gion III interphase as a whole of 2.197 GPa which is slightly less than the Case 2 value of 2.2039 GPa. 

Also of note in Figure 9, as the number of interphase subregions is increased, the difference between 
successive increases decreases. For example, the difference between the 2-Interphase case and the 3-Interphase 
case relative to the no CNT case is 1.66% at 80% carbon fiber volume fraction. At the same carbon fiber 
volume fraction, the difference between the 3-Interphase case and the 4-Interphase case is 0.34%. This 
demonstrates that there is a converging effect of adding increasing number of interphases which improves 
the accuracy with which the graded interphase region is modeled. As such, the comparison between the most 
accurate case, the 4-Interphase case, and the No CNT case is provided in Figure 10 to clearly demonstrate 
the effectiveness of adding even a small amount of CNTs on the lamina properties of carbon fiber lamina. 

IV. Conclusion 

A multi-scale method has been applied to assess the effective properties of a CNT coated carbon fiber 
lamina. The method makes use of the molecular dynamics, equivalent continuum, and micromechanics 
modeling approaches. Molecular dynamics simulations in conjunction with the equivalent continuum method 
have been used to obtain effective properties for a molecular RVE consisting of CNTs in epoxy. The effective 
CNT-epoxy properties are then randomized using the Mori- Tanaka method. Finally, a multi-layer composite 
cylinders approach has been applied to determine the effective properties of the carbon fiber lamina using 
an RVE consisting of up to four regions, two of which are interphase regions between the carbon fiber and 
the epoxy matrix consisting of randomly oriented CNTs in epoxy. The CNT volume fraction within these 
interphase regions has been taken to be variable resulting in a graded interphase region which is modeled using 
the multi-layer composite cylinders method using increasing numbers of piecewise continuous subregions. 



Figure 9. Effective plane strain bulk modulus (^ 23 ) of 
CNT-sprayed carbon fiber lamina of varying interphase 
representation. Results are shown only for the high car- 
bon fiber volume fractions for clarity. The labels NoCNTs, 
CNTs Near C-Fiber, and CNTs Far from C-Fiber have been 
identified in Figure 7(b). The new labels CNTs Near and 
Far from C-Fiber denote Lamina RVE representations 
including both Regions II and III, with the number 
in parenthesis denoting the total number of interphase 
regions. The order of the curves in the figure from bot- 
tom to top is: No CNTs, CNTs Near C-fiber (1), CNTs Near 
and Far from C-fiber (4), CNTs Near and Far from C-fiber 
(3), CNTs Far from C-fiber (1), CNTs Near and Far from 
C-fiber (2). 
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Effective lamina properties for various carbon 
fiber volume fractions were obtained, and it was 
found that the axial Young’s modulus and Poisson’s 
ratio were only moderately affected by the presence 
of the CNT graded interphase region at all volume 
fractions. In contrast, the transverse Young’s mod- 
ulus, plane strain bulk modulus, and the axial and 
transverse shear moduli were found to be very sensi- 
tive to the inclusion of the graded interphase region, 
especially at the high carbon fiber volume fractions 
typically found in composite lamina. For a lamina 
with 60% carbon fiber volume fraction, the results 
indicate that the transverse lamina properties may 
be enhanced by as much as 8%. This projected im- 
pact is better put into perspective by considering the 
total composite volume fraction of CNTs which for 
this increase in transverse properties is only 0.36% 

CNTs. Thus it is observed that a very small volume 
fraction of CNTs (0.36%) applied to the matrix surrounding carbon fibers in a lamina can have a measurable 
impact (8%) on the effective elastic properties associated with the direction transverse to the carbon fiber 
axis (E 2 2 , k 23 , G 2 3 and Gw). 



Figure 10. Effective plane strain bulk modulus (^ 23 ) of 
4-Interphase representation of a CNT-sprayed carbon 
fiber lamina. Results are shown only for the common 
carbon fiber volume fractions in a given lamina. 
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